JEF Paper Coding (2024) - v4.4 NYSE only (For Annual Data - Size & BM)

Notes to this new file: (March 26, 2024)

• This file builds on the file “Xu_2023_v4.4_AllfirmsPorts_reproduced.Rmd” and the original NYSE file “Xu_2022_v4.0.Rmd” (from folder “R&R_JEF_Mar31_2022”).

• The purpose of this replication/updated file is to ensure that my initial codes and results in the main paper are free from errors.

• Currently, the treatment of negative earnings is removing all negative earnings.

• Now the whole universe is the NYSE market, rather than the whole CRSP universe.

• For calculating the \(\bar{ge}\), I now update the formula.

• Variable vwret in the dataset is the log return.

• Variable freq now can be used to get results with different frequencies. Meanwhile, the results in the reply file to the AE and the table in the main paper should be updated.

Data cleaning and Preparation

timestamp() 

## ##------ Wed Mar 27 09:36:01 2024 ------##

# 0. record datasets ----
## 0.1 initial value setup ----
freq = 1 # the frequency of the data <- 12 for monthly; 4 for quarterly; 1 for annually
start.ym = as.yearmon(1966) -1/12 # the starting time
end.ym = as.yearmon(2020) # the ending time 
negativevalue = 'remove' # what should we do with the negative values? 
if (negativevalue == 'remove') cat("Negative earnings are deleted.") 

## Negative earnings are deleted.

month_select <- function(freq = freq) { # return the months for differnt time frequency
  if (freq == 12) {
    return(c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November","December"))
  }
  if (freq == 4) {
    # return(c("March", "June", "September", "December"))
    return(months(seq(from = as.yearmon(1966) + 2/12, by = 1/4, length.out = 4)))
  }
  if (freq == 1) {
    # return("December") # > here is unfavorable month 
    return(months(as.yearmon(1966) + 6/12)) # > July observations 
  }
}
freq_name <- function(freq = freq) {
  if (freq == 12) {
    return("monthly")
  }
  if (freq == 4) {
    return("quarterly") 
  }
  if (freq == 1) {
    return("annual")
  }
}

## 0.2 preliminary data cleaning ---- 
cat(paste("Selected month:", paste(month_select(freq = freq), collapse = ", ")), sep = "")

## Selected month: July

# "Variables_DK_CRSP_monthly_196601_201912.csv"
# "PredictorDataDK_CRSPm.csv"
# "PredictorDataDK_monthly_all.csv" - DK replication 
port_market <- read.csv("NYSE_compv6.csv") %>% 
  as_tibble() %>% 
  mutate(month = as.yearmon(as.Date(date)), 
         vwret = rollsumr(vwret, 12/freq, NA) ) %>% # converted into the log returns. 
  filter(if (negativevalue == 'remove') E >= 0 else TRUE) %>% 
  mutate(E = ((E*(E >= 0) + 0.001*(E < 0))) * 12/freq ) %>% 
  filter(months(month) %in% month_select(freq = freq)) %>% 
  filter(month < end.ym) 

plot(y = port_market$E, x = port_market$month, type = "b"); abline(h = 0, col = 'red', lty = 2)

    #elif negative == 'truncate':
    #    df['E'] = np.where(df['E'] < 0, 0.001, df['E'])

write.csv(x = port_market, file = "market_Allfirms.csv")
# port_market_raw <- port_market %>% select(month, vwret, R, date, GM, GE, DP)

ports_all <- read.csv("NYSE_portv3.csv") %>% 
  # PredictorDataDK_port_CRSPm.csv
  # PredictorDataDK_szbm_CRSPm.csv
  as_tibble() 

ports <- unique(ports_all$port) # identifiers for portfolios
for (p in ports) {
  ports.dt <- ports_all %>%
    filter(port == p) %>% 
    select(-port) %>% 
    mutate(month = as.yearmon(as.Date(date, "%d/%m/%Y")) ) %>% 
    arrange(month) %>% 
    mutate(vwret = rollsumr(vwret, 12/freq, NA) ) %>% 
    filter(if (negativevalue == 'remove') E >= 0 else TRUE) %>% 
    mutate(E = ((E*(E >= 0) + 0.001*(E < 0))) * 12/freq ) %>% 
    filter(months(month) %in% month_select(freq = freq)) %>% 
    filter(month < end.ym) 
  
  write.csv(x = ports.dt, file = paste("port_", p, ".csv", sep = ""))
}

## 0.3 name portfolios and predictors ----
market.names <- list.files(pattern = "market_")[1]
data.names <- list.files(pattern = "^port_") # data for portfolios
  ## generate the name of portfolios 
  ## Define the two sets
  #set1 <- c("B", "S")
  #set2 <- c("H", "M", "L")
  # Use expand.grid() to generate all combinations
  #combinations <- expand.grid(set2, set1)[, c(2, 1)]
  ## Convert the result to a vector
  # ports <- do.call(paste0, combinations) ## ----updated March 7, 2024 
id.names <- c("Market", ports) # set plot names
ratio_names <- c("DP", "PE", "EY", "DY", "Payout") # potential predictors

## 0.4 risk-free rate ----- 
RF <- read.csv("Rfree_t30.csv") %>% # record the risk free rate
  as.tbl() %>% # as the average of the bid and ask.
  select(-X) %>%
  mutate(month = as.yearmon(month)) %>%
  filter(months(month) %in% month_select(freq = freq)) %>%
  filter(month < end.ym)

Figure 1 - Log Cumulative Index

Log cumulative realised portfolio return components for seven portfolios - the market portfolio and six size and book-to-market equity ratio sorted portfolios. All following figures demonstrate the annual realized price-earnings ratio growth (gm), earnings growth (ge), dividend-price (dp) and the portfolio return index (r) with the values in January 1966 as zero for all portfolios.

Log Cumulative Index

# TABLE-1. summary statistics ----
TABLE1.uni <- list() # the univariate statistics
TABLE1.cor <- list() # the correlation matrixs

month.dt <- seq(filter(select(port_market, month), month >= start.ym)$month[2],
                tail(port_market$month, 1), 1/freq)
# head(month.dt)
DP.df <- EY.df <- PE.df <- data.frame(month = month.dt)

## (1*) summary tables for Summary stats & Correlations ----
for (c in 1:7) {
  id <- c(market.names, data.names)[c]
  # print(id); print(id.names[c])
  
  ## 1. read the data ---- 
  data_nyse <- read.csv(id) %>%
    as_tibble() %>%
    mutate(month = as.yearmon(month)) %>% 
    filter(month >= start.ym) %>% # start from "Dec 1965" 
    select(month, r = vwret, P, E = E, D = D, pe_raw = PE) %>% 
    # , R, GM_raw = GM, GE_raw = GE, DP_raw = DP
    mutate(DP = D / P, # these are adjusted by the log transformation
           PE = P / E,
           EP = E / P,
           EY = E / lag(P), # earnings yield
           DY = D / lag(P), # dividend yield
           Payout = D / E) # payout ratios
  
  PE.df <- PE.df %>% 
    left_join(select(data_nyse, month, PE), by = 'month')
  EY.df <- EY.df %>% 
    left_join(select(data_nyse, month, EY), by = 'month')
  DP.df <- DP.df %>% 
    left_join(select(data_nyse, month, DP), by = 'month')
  
  ## 2. return decomposition ----
  data_decompose <- data_nyse %>% 
    mutate(r = r, # cts returns = log total returns 
           gm = log(PE) - lag(log(PE)), # multiple expansion rate
           ge = log(E) - lag(log(E)), # earnings growth rate
           dp = log(1 + DP/freq)) %>% # only 1/12 of the dividends -----updated March 7, 2024 
    na.omit()
  
  ## 3. summary-Stat ----
  ar1.coef <- function(x) {
    return(as.numeric(lm(x ~ lag(x))$coefficients[2]))
  } # return the function value of the coefficient for the AR(1) model
  
  comp_summary.dt <- data_decompose %>%
    select(gm, ge, dp, r) %>%
    # , R, GM_raw, GE_raw, DP_raw 
    describe() %>%
    mutate(mean = mean * 100,
           sd = sd * 100,
           median = median * 100,
           min = min * 100,
           max = max * 100) %>%
    select(Mean = mean, Median = median, SD = sd, Min = min, Max = max, Skew = skew, Kurt = kurtosis) %>%
    round(digits = 4)
  
  comp_summary.dt$"AR(1)" <- data_decompose %>%
    select(gm, ge, dp, r) %>%
    apply(2, ar1.coef) %>%
    round(digits = 4)
  
  ### Store the summary stat
  # print(paste("Data starts from ", first(data_decompose$month), " and ends in ", last(data_decompose$month), ".", sep = ""))
  TABLE1.uni[[id.names[c]]] <- comp_summary.dt
  
  ## 4. correlations ----
  comp_cor <- data_decompose %>% select(gm, ge, dp, r) %>% cor()
  TABLE1.cor[[id.names[c]]] <- comp_cor
  
  # Figure-1. cumulative realised return components ---- 
  cat('\n')
  cat('#### ', id.names[c], '   \n')
  cat('\n')
  
  # jpeg(filename = paste("Figure1_", id.names[c], ".jpeg", sep = ""), width = 550, height = 350)
  par(mar = c(2, 4, 2, 1)) 
  cum_components.ts <- data.frame(month = month.dt) %>% 
    left_join(data_decompose, by = 'month') %>%
    select(r, gm, ge, dp) %>% 
    apply(2, FUN = function(x) cumsum(ifelse(is.na(x), 0, x) ) + x*0 ) %>% 
    ts(start = month.dt[1], frequency = freq)
  plot.ts(cum_components.ts, plot.type = "single", lty = 1:4, main = id.names[c], cex.main = 1, 
          xlab = NULL, ylab = "Cumulative Return and Components Indices")
  legend("topleft",
         legend = c("Total return", "Price earnings growth",
                    "Earnings growth", "Dividend price"),
         lty = 1:4,
         cex = 1.0) # text size
  # dev.off()
  par(mar = c(5, 4, 4, 2) + 0.1) 
  cat('\n')

} 

Market

BH

BL

BM

SH

SL

SM

write.csv(TABLE1.uni, file = "table_1.uni.csv")
write.csv(TABLE1.cor, file = "table_1.cor.csv")

.

Table 1 - Summary statistics of returns components

The correlations between gm and ge might be a bit too high comparing to Ferreira and Santa-Clara (2011). Need to check the code again.

‘kable’ for Table Creation

## summary data for table1
rowname <- rep(colnames(TABLE1.cor$Market), length(names(TABLE1.uni)))
portname = rep(names(TABLE1.uni), each = length(colnames(TABLE1.cor$Market)))

table1 <- do.call(rbind.data.frame, TABLE1.uni) %>% 
  cbind.data.frame(do.call(rbind.data.frame, TABLE1.cor)) %>%
  round(digits = 2) %>%
  cbind.data.frame(rowname, portname)

## give the table 1 outputs
gt(data = table1, rowname_col = "rowname", groupname_col = "portname") %>%
  tab_header(title = "Table 1 - Summary statistics of returns components",
             subtitle = paste(freq_name(freq = freq), " data starts from ", first(data_decompose$month), " and ends in ", last(data_decompose$month), ".", sep = "")) %>%
  tab_spanner(label = "Panel A: univariate statistics",
              columns = vars(Mean, Median, SD, Min, Max, Skew, Kurt, "AR(1)")) %>%
  tab_spanner(label = "Panel B: Correlations",
              columns = vars(gm, ge, dp, r)) %>%
  tab_source_note(source_note = html("Note: Panel A in this table presents mean, median, standard deviation (SD), minimum, maximum, skewness (Skew), kurtosis (kurt) and first-order autocorrelation coefficient of the realised components of stock market returns and six size and book-to-market equity ratio sorted portfolios. These univariate statistics for each portfolios are presented separately. gm is the continuously compounded growth rate in the price-earnings ratio. ge is the continuously compounded growth rate in earnings. dp is the log of one plus the dividend-price ratio. r is the portfolio returns. Panel B in this table reports correlation matrices for all seven portfolios. The sample period starts from Feburary 1966 and ends in December 2019."))

Table 1 - Summary statistics of returns components
annual data starts from Jul 1967 and ends in Jul 2019.
	Panel A: univariate statistics								Panel B: Correlations
	Mean	Median	SD	Min	Max	Skew	Kurt	AR(1)	gm	ge	dp	r
Market
gm	-0.09	0.01	37.78	-167.12	138.88	-0.67	8.06	-0.27	1.00	-0.93	-0.01	0.15
ge	7.67	10.78	38.85	-163.23	181.33	0.00	11.82	-0.20	-0.93	1.00	-0.05	0.23
dp	3.31	3.05	1.10	1.69	6.23	0.48	-0.71	0.87	-0.01	-0.05	1.00	-0.08
r	10.80	13.44	14.59	-27.07	49.98	-0.21	0.56	-0.08	0.15	0.23	-0.08	1.00
BH
gm	0.69	-12.25	96.14	-261.62	374.86	1.00	5.25	-0.36	1.00	-0.99	-0.06	0.38
ge	8.81	14.21	92.91	-380.08	244.50	-1.40	5.92	-0.31	-0.99	1.00	0.02	-0.24
dp	4.47	4.24	1.87	2.00	9.63	0.59	-0.41	0.75	-0.06	0.02	1.00	-0.12
r	12.57	13.23	14.97	-31.20	47.15	-0.51	0.99	-0.13	0.38	-0.24	-0.12	1.00
BL
gm	0.39	1.39	19.02	-59.92	47.57	-0.51	1.09	-0.19	1.00	-0.58	-0.22	0.71
ge	7.29	7.99	13.97	-42.26	51.52	-0.52	3.23	-0.17	-0.58	1.00	-0.04	0.15
dp	2.25	2.13	0.65	1.20	3.90	0.71	0.00	0.79	-0.22	-0.04	1.00	-0.22
r	10.35	11.82	15.64	-33.62	47.12	-0.36	0.26	-0.07	0.71	0.15	-0.22	1.00
BM
gm	0.40	2.00	34.13	-116.69	78.71	-0.66	1.78	-0.18	1.00	-0.90	-0.11	0.17
ge	7.60	9.84	34.72	-91.93	128.69	0.18	3.02	-0.09	-0.90	1.00	0.01	0.26
dp	3.59	3.06	1.39	1.75	7.40	0.88	-0.03	0.85	-0.11	0.01	1.00	-0.13
r	11.21	13.38	14.85	-20.89	50.20	-0.05	0.40	-0.09	0.17	0.26	-0.13	1.00
SH
gm	6.87	10.69	74.57	-325.56	162.15	-1.48	6.11	-0.01	1.00	-0.97	0.05	-0.03
ge	7.14	5.45	77.69	-163.73	360.37	1.53	6.88	0.04	-0.97	1.00	-0.12	0.25
dp	5.20	5.05	2.41	1.89	16.57	2.48	8.73	0.31	0.05	-0.12	1.00	-0.25
r	15.52	16.43	17.17	-15.57	66.76	0.15	-0.12	-0.09	-0.03	0.25	-0.25	1.00
SL
gm	-0.03	-5.10	91.77	-315.86	302.11	-0.68	4.68	-0.45	1.00	-0.97	0.21	0.23
ge	12.35	13.18	89.46	-301.10	325.89	0.62	5.49	-0.40	-0.97	1.00	-0.27	0.01
dp	2.79	2.18	2.46	0.62	14.12	3.04	10.07	0.43	0.21	-0.27	1.00	-0.18
r	12.41	13.72	20.63	-34.99	75.42	0.11	0.59	-0.25	0.23	0.01	-0.18	1.00
SM
gm	-0.42	-3.81	74.71	-294.34	298.48	0.24	7.12	-0.34	1.00	-0.97	-0.16	0.14
ge	13.02	15.30	74.71	-279.78	322.75	0.10	7.63	-0.35	-0.97	1.00	0.13	0.10
dp	3.69	3.31	1.74	1.51	9.29	1.35	1.71	0.33	-0.16	0.13	1.00	-0.05
r	14.50	17.33	17.79	-26.27	72.41	0.11	0.84	-0.20	0.14	0.10	-0.05	1.00
Note: Panel A in this table presents mean, median, standard deviation (SD), minimum, maximum, skewness (Skew), kurtosis (kurt) and first-order autocorrelation coefficient of the realised components of stock market returns and six size and book-to-market equity ratio sorted portfolios. These univariate statistics for each portfolios are presented separately. gm is the continuously compounded growth rate in the price-earnings ratio. ge is the continuously compounded growth rate in earnings. dp is the log of one plus the dividend-price ratio. r is the portfolio returns. Panel B in this table reports correlation matrices for all seven portfolios. The sample period starts from Feburary 1966 and ends in December 2019.

Figure 3 - Cumulative OOS R-sqaure Difference and Cumulative SSE Difference

Note:

• “SOP_c” here is for the constrained SOP predictions.

Cumulative OOS R-sqaure Difference and Cumulative SSE Difference

# TABLE-2. different tables OOS R-square ----
TABLE2 <- list()
table2.df <- data.frame()

## (3*) repeat for each portfolio ----
c <- 0
for (id in c(market.names, data.names)) {
  c <- c + 1
  # print(id); print(id.names[c])
  cat('\n')
  cat('#### ', id.names[c], '   \n')
  
  ## 1. read the data ----
  data_nyse <- read.csv(id) %>% 
    as.tbl() %>%
    mutate(month = as.yearmon(month)) %>%
    filter(month >= start.ym) %>% # start from "Dec 1965"
    select(month, r = vwret, P, E = E, D = D) %>%
    mutate(DP = D / P, # construct predictors
           PE = P / E,
           EP = E / P,
           EY = E / lag(P), # earnings yield
           DY = D / lag(P), # dividend yield
           Payout = D / E) # payout ratios
  
  ## 2. return decomposition ----
  k = freq * 20 # set a 20-year rolling window (total k periods.)
  
  data_decompose <- data_nyse %>% # also try PD ratio replacing PE.
    mutate(r = r, # cts returns (has already being the log return in row 95)
           gm = log(PE) - lag(log(PE)), # multiple expansion rate
           ge = log(E) - lag(log(E)), # earnings growth rate
           # mu_ge0 = (log(E) - lag(log(E), k)) / k,
           dp = log(1 + DP/freq)) %>% # see note 1.
    # only 1/12 of the annualised dividend is included. > see note 1.
    na.omit() %>% 
    left_join(select(data_nyse, month, Ek = E) %>% mutate(month = month + k/freq), by = 'month' ) %>% 
    mutate(mu_ge0 = (log(E) - log(Ek)) / k ) 
  
  ## 3. SOP predictions ---- 
  data_pred <- data_decompose %>%
    select(month, r, gm, ge, dp, mu_ge0) %>%
    mutate(mu_gm = 0,
           mu_ge1 = rollmeanr(ge, k, fill = NA_real_), # rolling mean 
           mu_ge2 = c(rep(NA_real_, (k-1)), cummean(ge)[k:length(ge)]), # recursive mean
           mu_ge3 = cummean(ge), # recursive mean from the beginning 
           a_DK1 = rollmeanr(r - dp, k, fill = NA_real_), # methods Eq (14/15) by DK 
           a_DK2 = cummean(r - dp), # methods Eq (14/15) by DK 
           mu_dp = dp,
           mu_sop = mu_gm + mu_ge0 + mu_dp) %>% # the predictor > see note 2.
    mutate(sop_simple = lag(mu_sop), # conditional predictions
           hist_mean = lag(cummean(r)) ) # historical mean predictions
  
  ### 3.2 constrained SOP predictions ----- 
    constrained_sop2 <- constrained_sop <- rep(NA_real_, nrow(data_decompose)) 
    for (t in (k+2):nrow(data_decompose)) {
      ## based on the Eq (14/15) in F-SC (2011) 
      x.IS <- data_decompose$dp[1:(t-2)] 
      y.IS <- data_decompose$r[2:(t-1)] 
      intercept <- data_pred$mu_ge1[t-1] #?
      y <- y.IS - intercept 
      data.IS <- data.frame(y, x.IS)
      
      model_constrained <- lm(y ~ x.IS - 1 + offset(1*x.IS), data = data.IS)
      model.IS <- glm(y ~ x.IS - 1, data = data.IS, family = gaussian(link = "identity"), offset = -1*x.IS)
      # mm <- glm.fit(x = x.IS, y = y, offset = -1*x.IS, intercept = F); mm$coefficients
      
      x.new <- data_decompose$dp[t-1] 
      y.pred <- predict.glm(model.IS, newdata = data.frame(x.IS = x.new), type = "response")
      # coef_constrained <- coef(model.IS) - 1 # the true coefficient estimate 
      y.pred_manual <- coef(model.IS) * x.new 
      constrained_sop[t] <- y.pred + intercept # store the prediction
      constrained_sop2[t] <- y.pred_manual + intercept # store the prediction
    }
    data_pred$sop_constrained <- constrained_sop
    data_pred$sop_constrained2 <- constrained_sop2
  
  ## 4. OOS R2 and MSE-F ----
  ### build the function for the bootstrap
  Boot_MSE.F <- function(data, actual, cond, uncond, x, critical.values = TRUE, boot.times = 10000) {
    ## clean and reorganise the data
    ports.dt <- data %>%
      select(month, actual = actual, cond = cond, uncond = uncond, x = x)
    
    ports.dt_narm <- ports.dt %>%
      na.omit() %>%
      mutate(error_A = actual - cond,
             error_N = actual - uncond,
             Cuml_MSE_A = cummean(error_A ^ 2),
             Cuml_MSE_N = cummean(error_N ^ 2),
             Cuml_OOS.rsquared = 1 - Cuml_MSE_A / Cuml_MSE_N,
             Cum_SSE = cumsum(error_N ^2 - error_A ^ 2) ) %>% 
      right_join(data.frame(month = seq(na.omit(ports.dt)$month[1], end.ym - 1/12, 1/freq) ), by = 'month') %>% 
      arrange(desc(-month))
    # ports.dt_narm %>% filter(month > 2009 & month < 2010.5)
    ### for the full (out-of-)sample
    MSE_A <- mean(ports.dt_narm$error_A ^ 2, na.rm = T)
    MSE_N <- mean(ports.dt_narm$error_N ^ 2, na.rm = T)
    OOS_r.squared <- 1 - MSE_A / MSE_N # out-of-sample R square
    MSE_F <- length(na.omit(ports.dt_narm)$month) * (MSE_N - MSE_A) / (MSE_A) 
    MAE_A <- mean(abs(ports.dt_narm$error_A), na.rm = T) %>% round(digits = 6) # Mean absolute error of the conditional model
    print(paste("OOS R Squared: ", round(OOS_r.squared, digits = 4), sep = ""))
    print(paste("MSE-F: ", round(MSE_F, digits = 4), sep = ""))
    
    Cuml_OOS.ts <- ts(ports.dt_narm$Cuml_OOS.rsquared, start = ports.dt_narm$month[1], frequency = freq)
    Cum_SSE.ts <- ts(ports.dt_narm$Cum_SSE, start = ports.dt_narm$month[1], frequency = freq)
    
    ## Bootstrapped MSE-F Stat
    if (critical.values == TRUE) {
      ## get the data series for x and y
      y0 <- ports.dt[is.na(ports.dt$cond), ]$actual
      y1 <- ports.dt[!is.na(ports.dt$cond), ]$actual
      x1 <- as.vector(na.omit(ports.dt$x))
      
      ## full sample estimation for the model
      alpha <- mean(y1)
      u1 <- y1 - alpha
      
      x.lm <- lm(x1 ~ lag(x1))
      mu <- as.numeric(coef(x.lm)[1])
      rho <- as.numeric(coef(x.lm)[2])
      u2 <- as.numeric(x.lm$residuals)
      
      u <- data.frame(index = 1:length(u1), u1, u2) # the dataset storing all the residual info
      
      ### bootstrapping pairs of error terms
      boot.critical <- c()
      for (i in 1:boot.times) { # the bootstrapped times can be modified
        index.boot <- sample(u$index, length(u1), replace = T)
        u.boot <- data.frame()
        for (j in index.boot) {
          u.boot <- rbind.data.frame(u.boot, u[j,])
        } # record the bootstrapped error terms
        
        ### reconstruct the simulated x and y
        y1.new <- alpha + u.boot$u1
        x0.new <- sample(x1, 1) # resample a value as the starting value of x
        x1.new <- c(mu + rho * x0.new + u.boot$u2[1])
        for (j in 2:length(y1.new)) {
          x1.new[j] <- mu + x1.new[j-1] * rho + u.boot$u2[j]
        } # simulate SOP values
        
        ### redo the rolling estimation
        y.boot <- c(y0, y1.new)
        x.boot <- c(rep(NA_real_, (length(y0) - 1)), x0.new, x1.new)
        
        data.dt <- data.frame(month = ports.dt$month, x.boot, y.boot) %>%
          as.tbl() %>%
          mutate(conditional = lag(x.boot), # convert to the SOP prediction
                 unconditional = lag(cummean(y.boot))) %>% # convert to the historical mean prediction
          na.omit()
        
        error_N.boot = data.dt$y.boot - data.dt$unconditional
        error_A.boot = data.dt$y.boot - data.dt$conditional
        MSE_N.boot = mean(error_N.boot ^ 2) 
        MSE_A.boot = mean(error_A.boot ^ 2) 
        OOS_r.squared.boot <- 1 - MSE_A.boot / MSE_N.boot %>% round(digits = 6) # out-of-sample R square
        
        ### MSE-F statistic
        MSE_F.boot <- length(error_N.boot) * (MSE_N.boot - MSE_A.boot) / (MSE_A.boot) %>% round(digits = 6)
        
        boot.critical[i] <- MSE_F.boot
        
        if (i %% (boot.times/10) == 0) {
          timestamp()
          print(paste("SOP", ": ", i %/% (boot.times/100), "%", sep = ""))
        }
      }
      ## store the results
      result <- cbind.data.frame(IS_r.squared = NA_real_,
                                 OOS_r.squared,
                                 MAE_A, # MAE of conditional model
                                 MSE_F,
                                 t(quantile(boot.critical, probs = c(0.90, 0.95, 0.99))),
                                 p.value = mean(boot.critical > MSE_F))
      # d_RMSE = round(d_RMSE, digits = 4),
      # DM_stat = round(DM_test.result$statistic, digits = 4),
      # DM_pval = round(DM_test.result$p.value, digits = 4))
    } else {
      result <- cbind.data.frame(IS_r.squared = NA_real_,
                                 OOS_r.squared,
                                 MAE_A, # MAE of conditional model
                                 MSE_F)
    }
    
    rownames(result) <- "SOP"
    output <- list(result = result, Cuml_OOS.ts = Cuml_OOS.ts, Cum_SSE.ts = Cum_SSE.ts)
    return(output)
  }
  
  ### store the results for the SOP method
  cat("SOP: ")
  sop.result <- Boot_MSE.F(data = data_pred, actual = "r", cond = "sop_simple", uncond = "hist_mean", x = "mu_sop", critical.values = FALSE, boot.times = 3000)
  sop.result$result 
  cat('\n')
  
  cat("SOP_c: ")
  sop_constrained.result <- Boot_MSE.F(data = data_pred, actual = "r", cond = "sop_constrained", uncond = "hist_mean", x = "mu_sop", critical.values = FALSE, boot.times = 3000)
  sop_constrained.result$result
  cat('\n')
  
  # sop_constrained2.result <- Boot_MSE.F(data = data_pred, actual = "r", cond = "sop_constrained2", uncond = "hist_mean", x = "mu_sop", critical.values = FALSE, boot.times = 3000)
  # sop_constrained.result$result
  
  cat('\n')
  par(mfrow = c(2,1))
    par(mar = c(3,5,3,2))
  Cuml_OOS.sop.ts <- ts.intersect(sop.result$Cuml_OOS.ts * 100, sop_constrained.result$Cuml_OOS.ts * 100) 
  plot.ts(Cuml_OOS.sop.ts, 
          plot.type = 'single', lty = c(1,3), col = c(1,'blue'), 
          xlab = NULL, ylab = "as %", main = paste(id.names[c], ": Cumulative OOS R Squared Difference - SOP", sep = "")) 
  abline(h = c(0,1), lty = 2, col = 2)
  abline(v = c(2003,2007), lty = 2, col = 2)
  
    par(mar = c(5,5,3,2))
  Cuml_SSE.sop.ts <- ts.intersect(sop.result$Cum_SSE.ts, sop_constrained.result$Cum_SSE.ts) 
  plot.ts(Cuml_SSE.sop.ts, 
          plot.type = 'single', lty = c(1,3), col = c(1,'blue'), 
          ylab = NULL, cex.lab = 0.5, 
          xlab = "An increase in a line indicates better performance of the conditional model\n *The blue dotted line is for the constrained SOP",
          main = paste(id.names[c], ": Cumulative SSE Difference - SOP", sep = "")) 
  abline(h = 0, lty = 2, col = 2)
  abline(v = c(2003,2008), lty = 2, col = 2)
  par(mfrow = c(1,1))
  
  cat('\n') 
  
  ### store all cumulative OOS R2 difference values
  Cuml_all.ts <- sop.result$Cuml_OOS.ts * 100
  Csse_all.ts <- sop.result$Cum_SSE.ts
  
  ## 5. univariate predictive regressions ----
  table2.uni_predictors <- data.frame()
  cat('\n')
  
  for (predictor in ratio_names) {
    ## construct conditional & unconditional predictions
    data_univariate <- data_decompose %>%
      select(month, r, predictor) %>%
      mutate(hist_mean = lag(cummean(r)),
             x = lag(get(predictor)) %>% log) ## convert to log predictors
      # data_univariate[[predictor]] <- log(data_univariate[[predictor]])
      # data_univariate[["x"]] <- log(data_univariate[["x"]])
    
    ## IS R2
    lm.IS <- lm(r ~ x, data = data_univariate)
    IS_r.squared <- summary(lm.IS)$r.squared # in-sample r squared
    IS_pval <- summary(lm.IS)$coefficients[2,4] # the p-value from F-statistic
    
    ## OOS recursive window predictions
    k <- k # the starting in-sample data
    con_pred = rep(NA_real_, nrow(data_univariate))
    for (t in (k+2):nrow(data_univariate)) {
      x.IS <- data_univariate$x[2:(t-1)]
      y.IS <- data_univariate$r[2:(t-1)]
      reg.IS <- lm(y.IS ~ x.IS)
      
      x.new <- data_univariate$x[t]
      y.pred <- predict(reg.IS, newdata = data.frame(x.IS = x.new))
      con_pred[t] <- y.pred # store the prediction
    }
    data_univariate$con_pred <- con_pred
    data_univariate
    
    ## Stat and Bootstrap
    data = data_univariate
    actual = "r"
    cond = "con_pred"
    uncond = "hist_mean"
    critical.values = FALSE # decide whether the bootstrapped critical value is calculated > Note 3.
    boot.times = 1000 * 10
    
    { 
      cat('\n')
      cat(paste(predictor, ": ")) 
      ## get OOS R2 & MSE-F
      ports.dt <- data %>%
        select(month, actual = actual, cond = cond, uncond = uncond, predictor) %>%
        rename(x = predictor)
      
      ports.dt_narm <- ports.dt %>%
      na.omit() %>%
      mutate(error_A = actual - cond,
             error_N = actual - uncond,
             Cuml_MSE_A = cummean(error_A ^ 2),
             Cuml_MSE_N = cummean(error_N ^ 2),
             Cuml_OOS.rsquared = 1 - Cuml_MSE_A / Cuml_MSE_N,
             Cum_SSE = cumsum(error_N ^2 - error_A ^ 2) ) %>% 
      right_join(data.frame(month = seq(na.omit(ports.dt)$month[1], end.ym - 1/12, 1/freq) ), by = 'month') %>% 
      arrange(desc(-month))
      
      ### for the full (out-of-)sample
      MSE_A <- mean(ports.dt_narm$error_A ^ 2, na.rm = T)
      MSE_N <- mean(ports.dt_narm$error_N ^ 2, na.rm = T)
      OOS_r.squared <- 1 - MSE_A / MSE_N # out-of-sample R square
      MSE_F <- length(na.omit(ports.dt_narm)$month) * (MSE_N - MSE_A) / (MSE_A) 
      MAE_A <- mean(abs(ports.dt_narm$error_A), na.rm = T) %>% round(digits = 6) # Mean absolute error of the conditional model
      print(paste("IS R Squared: ", round(IS_r.squared, digits = 4), sep = ""))
      print(paste("OOS R Squared: ", round(OOS_r.squared, digits = 4), sep = ""))
      print(paste("MSE-F: ", round(MSE_F, digits = 4), sep = ""))
      
      ### combine the cumulative OOS R2 difference of other predictors & cum SSE
      Cuml_pred.ts <- ts(ports.dt_narm$Cuml_OOS.rsquared * 100, start = ports.dt_narm$month[1], frequency = freq)
      Cuml_all.ts <- ts.union(Cuml_all.ts, Cuml_pred.ts)
      Cum_pred.sse <- ts(ports.dt_narm$Cum_SSE, start = ports.dt_narm$month[1], frequency = freq)
      Csse_all.ts <- ts.union(Csse_all.ts, Cum_pred.sse)
      
      ## prepare the bootstrap
      if (critical.values == TRUE) {
        y0 <- ports.dt$actual[1]
        y1 <- ports.dt$actual[-1]
        x1 <- ports.dt$x
        
        ## full sample estimation for the model
        alpha <- mean(y1)
        u1 <- y1 - alpha
        
        x.lm <- lm(x1 ~ lag(x1))
        mu <- as.numeric(coef(x.lm)[1])
        rho <- as.numeric(coef(x.lm)[2])
        u2 <- as.numeric(x.lm$residuals)
        
        u <- data.frame(index = 1:length(u1), u1, u2) # the dataset storing all the residual info
        
        ### bootstrap pairs of error terms
        boot.critical <- c()
        for (i in 1:boot.times) { # the bootstrapped times can be modified
          index.boot <- sample(u$index, length(u1), replace = T)
          u.boot <- data.frame()
          for (j in index.boot) {
            u.boot <- rbind.data.frame(u.boot, u[j,])
          } # record the bootstrapped error terms
          
          ### reconstruct the simulated x and y
          y1.new <- alpha + u.boot$u1
          x0.new <- sample(x1, 1) # resample a value as the starting value of x
          x1.new <- c(mu + rho * x0.new + u.boot$u2[1])
          for (j in 2:length(y1.new)) {
            x1.new[j] <- mu + x1.new[j-1] * rho + u.boot$u2[j]
          } # simulate predictors
          
          ### redo the rolling estimation
          y.boot <- c(y0, y1.new)
          x.boot <- c(x0.new, x1.new)
          
          data.dt <- as.tbl(data.frame(month = ports.dt$month, x.boot, y.boot, x = lag(x.boot)))
          con_pred = rep(NA_real_, nrow(data.dt))
          for (t in (k+2):nrow(data.dt)) {
            x.IS <- data.dt$x[2:(t-1)]
            y.IS <- data.dt$y.boot[2:(t-1)]
            reg.IS <- lm(y.IS ~ x.IS)
            
            x.new <- data.dt$x[t]
            y.pred <- predict(reg.IS, newdata = data.frame(x.IS = x.new))
            con_pred[t] <- y.pred
          }
          data.dt$conditional <- con_pred
          data.dt <- data.dt %>%
            mutate(unconditional = lag(cummean(y.boot))) %>%
            na.omit()
          
          error_N.boot = data.dt$y.boot - data.dt$unconditional
          error_A.boot = data.dt$y.boot - data.dt$conditional
          MSE_N.boot = mean(error_N.boot ^ 2)
          MSE_A.boot = mean(error_A.boot ^ 2)
          OOS_r.squared.boot <- 1 - MSE_A.boot / MSE_N.boot # out-of-sample R square
          
          ### MSE-F statistic
          MSE_F.boot <- length(error_N.boot) * (MSE_N.boot - MSE_A.boot) / (MSE_A.boot)
          
          boot.critical[i] <- MSE_F.boot
          
          if (i %% (boot.times/10) == 0) {
            timestamp()
            print(paste(predictor, ": ", i %/% (boot.times/100), "%", sep = ""))
          }
        }
        ## store the results
        result <- cbind.data.frame(IS_r.squared, 
                                   OOS_r.squared,
                                   MAE_A, # MAE of conditional model
                                   MSE_F,
                                   t(quantile(boot.critical, probs = c(0.90, 0.95, 0.99))),
                                   p.value = mean(boot.critical > MSE_F))
        # d_RMSE = round(d_RMSE, digits = 4),
        # DM_stat = round(DM_test.result$statistic, digits = 4),
        # DM_pval = round(DM_test.result$p.value, digits = 4))
      } else {
        result <- cbind.data.frame(IS_r.squared, 
                                   OOS_r.squared,
                                   MAE_A, # MAE of conditional model
                                   MSE_F)
      }
      rownames(result) <- predictor
      
      result 
      cat('\n')
    }
    
    ## store the results for all predictors
    table2.uni_predictors <- rbind.data.frame(table2.uni_predictors, result)
  }
  table2.uni_predictors
  
  ## 6. statistics summary ----
  table2 <- rbind.data.frame(table2.uni_predictors, "SOP" = sop.result$result, "SOP_c" = sop_constrained.result$result)
  if ("p.value" %in% colnames(table2)) {
    table2$star <- ifelse(table2$p.value <= 0.01, "***", ifelse(table2$p.value <= 0.05, "**", ifelse(table2$p.value <= 0.1, "*", "")))
  }
  # statistical significance from McCracken (2007)
  table2$McCracken <- ifelse(table2$MSE_F >= 3.838, "***", ifelse(table2$MSE_F >= 1.599, "**", ifelse(table2$MSE_F >= 0.685, "*", "")))
  
  table2
  TABLE2[[id.names[c]]] <- table2
  
  table2$rowname <- rownames(table2)
  table2$portname <- id.names[c]
  table2.df <- rbind.data.frame(table2.df, table2)
  
  ## 7. cumulative OOS R2 difference merged dataset ----
  colnames(Cuml_all.ts) <- c("SOP", ratio_names)
  cat('\n') 
  par(mfrow = c(3,3))
  for (method in colnames(Cuml_all.ts)) {
    par(mar = c(2, 2.5, 2, 0.5))
    plot.ts(window(Cuml_all.ts[, method], start = 1990), xlab = NULL, lty = 2, ylab = NULL, main = method, cex.main = 0.8, xaxt = "n", las = 2)
    axis(1, seq(1990, 2020, by = 10), cex.axis = 0.8)
  }
  
  colnames(Csse_all.ts) <- c("SOP", ratio_names)
  par(mfrow = c(3,3))
  for (method in colnames(Csse_all.ts)) {
    par(mar = c(2, 2.5, 2, 0.5))
    plot.ts(Csse_all.ts[, method], xlab = NULL, lty = 2, ylab = NULL, main = method, cex.main = 0.8, xaxt = "n", las = 2)
    axis(1, seq(1990, 2020, by = 10), cex.axis = 0.8)
  }
  par(mfrow = c(1,1), mar = c(5, 4, 4, 2) + 0.1)
  cat('\n')
}

Market

SOP: [1] “OOS R Squared: 0.0383” [1] “MSE-F: 1.3139”

SOP_c: [1] “OOS R Squared: 0.0409” [1] “MSE-F: 1.3636”

DP : [1] “IS R Squared: 0.1088” [1] “OOS R Squared: -0.2719” [1] “MSE-F: -6.8402”

PE : [1] “IS R Squared: 0.0471” [1] “OOS R Squared: -0.0711” [1] “MSE-F: -2.1253”

EY : [1] “IS R Squared: 0.0246” [1] “OOS R Squared: -0.0689” [1] “MSE-F: -2.0631”

DY : [1] “IS R Squared: 0.0705” [1] “OOS R Squared: -0.0424” [1] “MSE-F: -1.3011”

Payout : [1] “IS R Squared: 0.0048” [1] “OOS R Squared: -0.0747” [1] “MSE-F: -2.2242”

BH

SOP: [1] “OOS R Squared: 0.1958” [1] “MSE-F: 8.0368”

SOP_c: [1] “OOS R Squared: 0.1741” [1] “MSE-F: 6.7469”

DP : [1] “IS R Squared: 0.0954” [1] “OOS R Squared: -0.0421” [1] “MSE-F: -1.2916”

PE : [1] “IS R Squared: 0.1654” [1] “OOS R Squared: 0.0877” [1] “MSE-F: 3.0758”

EY : [1] “IS R Squared: 0.145” [1] “OOS R Squared: 0.0307” [1] “MSE-F: 1.0148”

DY : [1] “IS R Squared: 0.0684” [1] “OOS R Squared: -0.0293” [1] “MSE-F: -0.9097”

Payout : [1] “IS R Squared: 0.0867” [1] “OOS R Squared: -0.167” [1] “MSE-F: -4.5799”

BL

SOP: [1] “OOS R Squared: 0.0184” [1] “MSE-F: 0.6178”

SOP_c: [1] “OOS R Squared: 0.0207” [1] “MSE-F: 0.6748”

DP : [1] “IS R Squared: 0.1129” [1] “OOS R Squared: -0.0285” [1] “MSE-F: -0.8872”

PE : [1] “IS R Squared: 0.0689” [1] “OOS R Squared: -0.0262” [1] “MSE-F: -0.8161”

EY : [1] “IS R Squared: 0.0414” [1] “OOS R Squared: 0.0084” [1] “MSE-F: 0.2726”

DY : [1] “IS R Squared: 0.0766” [1] “OOS R Squared: 0.0501” [1] “MSE-F: 1.6893”

Payout : [1] “IS R Squared: 0.0094” [1] “OOS R Squared: -0.1397” [1] “MSE-F: -3.9218”

BM

SOP: [1] “OOS R Squared: -0.0038” [1] “MSE-F: -0.1262”

SOP_c: [1] “OOS R Squared: 0.023” [1] “MSE-F: 0.753”

DP : [1] “IS R Squared: 0.0393” [1] “OOS R Squared: -0.1943” [1] “MSE-F: -5.2055”

PE : [1] “IS R Squared: 0.0574” [1] “OOS R Squared: -0.0739” [1] “MSE-F: -2.2033”

EY : [1] “IS R Squared: 0.034” [1] “OOS R Squared: -0.0494” [1] “MSE-F: -1.5053”

DY : [1] “IS R Squared: 0.0226” [1] “OOS R Squared: -0.085” [1] “MSE-F: -2.5063”

Payout : [1] “IS R Squared: 0.0091” [1] “OOS R Squared: -0.0441” [1] “MSE-F: -1.3527”

SH

SOP: [1] “OOS R Squared: -0.1298” [1] “MSE-F: -3.7906”

SOP_c: [1] “OOS R Squared: -0.0974” [1] “MSE-F: -2.8391”

DP : [1] “IS R Squared: 0.0369” [1] “OOS R Squared: -0.1354” [1] “MSE-F: -3.8158”

PE : [1] “IS R Squared: 0.0064” [1] “OOS R Squared: -0.0361” [1] “MSE-F: -1.114”

EY : [1] “IS R Squared: 0.0035” [1] “OOS R Squared: -0.0379” [1] “MSE-F: -1.1699”

DY : [1] “IS R Squared: 0.0179” [1] “OOS R Squared: -0.0458” [1] “MSE-F: -1.4023”

Payout : [1] “IS R Squared: 3e-04” [1] “OOS R Squared: -0.0456” [1] “MSE-F: -1.3953”

SL

SOP: [1] “OOS R Squared: -0.1613” [1] “MSE-F: -4.5827”

SOP_c: [1] “OOS R Squared: -0.427” [1] “MSE-F: -9.5755”

DP : [1] “IS R Squared: 0.04” [1] “OOS R Squared: -0.8677” [1] “MSE-F: -14.8668”

PE : [1] “IS R Squared: 0.0105” [1] “OOS R Squared: -0.3825” [1] “MSE-F: -8.8529”

EY : [1] “IS R Squared: 0.0014” [1] “OOS R Squared: -0.1484” [1] “MSE-F: -4.1347”

DY : [1] “IS R Squared: 0.0109” [1] “OOS R Squared: -0.1388” [1] “MSE-F: -3.8998”

Payout : [1] “IS R Squared: 6e-04” [1] “OOS R Squared: -0.1089” [1] “MSE-F: -3.1417”

SM

SOP: [1] “OOS R Squared: -0.3014” [1] “MSE-F: -7.6418”

SOP_c: [1] “OOS R Squared: -0.1362” [1] “MSE-F: -3.8355”

DP : [1] “IS R Squared: 0.0165” [1] “OOS R Squared: -0.7863” [1] “MSE-F: -14.0857”

PE : [1] “IS R Squared: 0.0017” [1] “OOS R Squared: -0.0038” [1] “MSE-F: -0.1221”

EY : [1] “IS R Squared: 3e-04” [1] “OOS R Squared: 0.0039” [1] “MSE-F: 0.1246”

DY : [1] “IS R Squared: 0.001” [1] “OOS R Squared: -0.2987” [1] “MSE-F: -7.3597”

Payout : [1] “IS R Squared: 0.0014” [1] “OOS R Squared: -0.0116” [1] “MSE-F: -0.3684”

Table 2 - Forecasts of portfolio returns

This table demonstrates the in-sample and out-of-sample R-squares for the market and six size and book-to-market equity ratio sorted portfolios from predictive regressions and the Sum-of-the-Parts method. IS R-squares are estimated using the whole sample period and the OOS R-squares are calculated compare the forecast error of the model against the historical mean model. The full sample period starts from Feb 1966 to December 2019 and the IS period is set to be 20 years with forecsats beginning in Feb 1986. The MSE-F statistics are calculated to test the hypothesis \(H_0: \text{out-of-sample R-squares} = 0\) vs \(H_1: \text{out-of-sample R-squares} \neq 0\).

Predictors here are all in log terms.

gt(table2.df, rowname_col = "rowname", groupname_col = "portname") %>%
  tab_header(title = "Table 2 - Forecasts of portfolio returns",
             subtitle = paste("OOS ", freq_name(freq = freq), " data starts from ", first(data_decompose$month) + k/freq, " and ends in ", last(data_decompose$month), ".", sep = "")) %>%
  fmt_number(columns = 1:4, decimals = 6, suffixing = TRUE)

Table 2 - Forecasts of portfolio returns
OOS annual data starts from Jul 1987 and ends in Jul 2019.
	IS_r.squared	OOS_r.squared	MAE_A	MSE_F	McCracken
Market
DP	0.108799	−0.271872	0.102958	−6.840235
PE	0.047085	−0.071140	0.099337	−2.125281
EY	0.024637	−0.068914	0.097364	−2.063075
DY	0.070528	−0.042382	0.094677	−1.301072
Payout	0.004826	−0.074698	0.089981	−2.224182
SOP	NA	0.038290	0.090928	1.313865	*
SOP_c	NA	0.040872	0.085645	1.363630	*
BH
DP	0.095371	−0.042061	0.129317	−1.291612
PE	0.165361	0.087690	0.114749	3.075806	**
EY	0.145035	0.030738	0.115630	1.014811	*
DY	0.068358	−0.029259	0.121743	−0.909662
Payout	0.086737	−0.167029	0.127451	−4.579933
SOP	NA	0.195844	0.119093	8.036833	***
SOP_c	NA	0.174127	0.113120	6.746885	***
BL
DP	0.112940	−0.028514	0.094286	−0.887166
PE	0.068886	−0.026172	0.099143	−0.816143
EY	0.041415	0.008447	0.097703	0.272610
DY	0.076638	0.050143	0.091223	1.689288	**
Payout	0.009383	−0.139676	0.093610	−3.921832
SOP	NA	0.018376	0.096416	0.617768
SOP_c	NA	0.020651	0.088651	0.674754
BM
DP	0.039300	−0.194273	0.101348	−5.205453
PE	0.057358	−0.073945	0.098738	−2.203311
EY	0.033981	−0.049363	0.097445	−1.505304
DY	0.022641	−0.084976	0.094768	−2.506263
Payout	0.009120	−0.044137	0.096308	−1.352683
SOP	NA	−0.003840	0.095617	−0.126245
SOP_c	NA	0.022990	0.088615	0.752978	*
SH
DP	0.036901	−0.135389	0.128434	−3.815835
PE	0.006375	−0.036068	0.131557	−1.114003
EY	0.003549	−0.037946	0.129854	−1.169891
DY	0.017897	−0.045830	0.124686	−1.402285
Payout	0.000332	−0.045589	0.125232	−1.395252
SOP	NA	−0.129775	0.134336	−3.790634
SOP_c	NA	−0.097359	0.132067	−2.839072
SL
DP	0.039979	−0.867716	0.128204	−14.866777
PE	0.010526	−0.382461	0.127013	−8.852878
EY	0.001424	−0.148383	0.119221	−4.134742
DY	0.010885	−0.138782	0.114175	−3.899790
Payout	0.000633	−0.108865	0.116086	−3.141661
SOP	NA	−0.161264	0.119452	−4.582687
SOP_c	NA	−0.427009	0.120224	−9.575467
SM
DP	0.016517	−0.786282	0.133490	−14.085696
PE	0.001720	−0.003830	0.110039	−0.122107
EY	0.000308	0.003879	0.108690	0.124596
DY	0.001030	−0.298685	0.120964	−7.359693
Payout	0.001375	−0.011645	0.109071	−0.368355
SOP	NA	−0.301357	0.127595	−7.641846
SOP_c	NA	−0.136183	0.115814	−3.835511

Figure 4 - Annual return predictions

Here I only present the annual predictions of the historical mean model, the SOP method and the predictive regressions based on the dividend-price ratio and the earnings-price ratio.

Return Predictions

# TABLE-3 optimal weights and CEGs ----
all_predictions <- list() ## store all the prediction values
# kpss.ts <- kpss.diff <- as.data.frame(matrix(nrow = 8, ncol = 7), col.names = c(market.names, data.names))
c <- 0

for (id in c(market.names, data.names)) {
  c <- c + 1
  # print(id); print(id.names[c])
  cat('\n')
  cat('#### ', id.names[c], '   \n')
  
  ## construct SOP predictions
  data_nyse <- read.csv(id) %>%
    as.tbl() %>%
    mutate(month = as.yearmon(month)) %>%
    filter(month >= start.ym) %>% # start from "Jan 1966"
    select(month, r = vwret, P, E = E, D = D) %>%
    mutate(DP = D / P,
           PE = P / E,
           EP = E / P,
           EY = E / lag(P), # earnings yield
           DY = D / lag(P), # dividend yield
           Payout = D / E) # payout ratios
  
  ## 2. return decomposition ----
  k = freq * 20 # set a 20-year rolling window 
  
  data_decompose <- data_nyse %>% # also try PD ratio replacing PE.
    mutate(r = r, # cts returns (has already being the log return in row 95)
           gm = log(PE) - lag(log(PE)), # multiple expansion rate
           ge = log(E) - lag(log(E)), # earnings growth rate
           # mu_ge0 = (log(E) - lag(log(E), k)) / k,
           dp = log(1 + DP/freq)) %>% # see note 1.
    # only 1/12 of the annualised dividend is included. > see note 1.
    na.omit() %>% 
    left_join(select(data_nyse, month, Ek = E) %>% mutate(month = month + k/freq), by = 'month' ) %>% 
    mutate(mu_ge0 = (log(E) - log(Ek)) / k ) 
  
  ## 3. SOP predictions ---- 
  data_pred <- data_decompose %>%
    select(month, r, gm, ge, dp, mu_ge0) %>%
    mutate(mu_gm = 0,
           mu_ge1 = rollmeanr(ge, k, fill = NA_real_), # rolling mean 
           mu_ge2 = c(rep(NA_real_, (k-1)), cummean(ge)[k:length(ge)]), # recursive mean
           mu_ge3 = cummean(ge), # recursive mean from the beginning 
           a_DK1 = rollmeanr(r - dp, k, fill = NA_real_), # methods Eq (14/15) by DK 
           a_DK2 = cummean(r - dp), # methods Eq (14/15) by DK 
           mu_dp = dp,
           mu_sop = mu_gm + mu_ge0 + mu_dp) %>% # the predictor > see note 2.
    mutate(sop_simple = lag(mu_sop), # conditional predictions
           hist_mean = lag(cummean(r)) ) # historical mean predictions
  
  all_predictions.dt <- data_pred %>% 
    right_join(data.frame(month = seq(data_pred$month[1], end.ym - 1/12, 1/freq) ), by = 'month') %>% 
    arrange(desc(-month)) %>% 
    select(month, r, hist_mean, sop_simple) # store the prediction results
  
  ## construct univariate predictions
  for (predictor in ratio_names) {
    ## construct conditional & unconditional predictions
    data_univariate <- data_decompose %>%
      select(month, r, predictor) %>%
      mutate(x = lag(get(predictor)) %>% log) %>%
      select(-predictor)
    
    ## IS R2
    lm.IS <- lm(r ~ x, data = data_univariate)
    IS_r.squared <- summary(lm.IS)$r.squared # in-sample r squared
    
    ## OOS recursive window predictions
    k <- k # the starting in-sample data
    con_pred = rep(NA_real_, nrow(data_univariate))
    for (t in (k+2):nrow(data_univariate)) {
      x.IS <- data_univariate$x[2:(t-1)]
      y.IS <- data_univariate$r[2:(t-1)]
      reg.IS <- lm(y.IS ~ x.IS)
      
      x.new <- data_univariate$x[t]
      y.pred <- predict(reg.IS, newdata = data.frame(x.IS = x.new))
      con_pred[t] <- y.pred # store the prediction
    }
    data_univariate[[predictor]] <- con_pred
    data_univariate # get the univariate predictions
    
    all_predictions.dt <- all_predictions.dt %>% # combine into all other estimations
      left_join(select(data_univariate, month, predictor), by = "month")
  }
  
  ## store all the predictions generated
  all_predictions[[id.names[c]]] <- na.omit(all_predictions.dt)
  names(all_predictions.dt)
  
 ## show the performance of variance predictions
  cat('\n')
  
  p <- c("hist_mean", "sop_simple", "DP", "PE")
  # all_predictions.ts <- na.omit(all_predictions.dt)
  all_predictions.ts <- ts(all_predictions.dt[, -(1:2)], start = all_predictions.dt$month[1], frequency = freq) %>% 
    window(start = as.numeric(start.ym + k/freq) )
  plot.ts(all_predictions.ts[, p], plot.type = "single",  col = 1:length(p), xlab = NULL, ylab = paste(str_to_title(freq_name(freq = freq)), " Returns", sep = ""))
  title(main = paste(id.names[c], ": ", freq_name(freq = freq), " return predictions", sep = ""))
  legend(xpd = T, inset = -0.35, "bottom", legend = p, col = 1:length(p), lty = 1, horiz = T, cex = 0.8, bty = "n", lwd = 2)
  cat('\n')
}

Market

BH

BL

BM

SH

SL

SM

Figure 5 - Trading Performance (with no trading restrictions)

• Markowitz optimal weight on the risky portfolio using the SOP method and the historical mean model.
• Markowitz optimal weight on the risky portfolio using all models.
• Out-of-sample performance of the trading strategies - the SOP method, the historical mean model and the buy-and-hold strategy - with no trading restrictions.

Trading Performance

## 8. calculate the CEGs ----
TABLE3 <- list()
TABLE4 <- list()

table3.df <- data.frame()
table4.df <- data.frame()

##
risk.coef <- 3 # the risk-aversion coefficient
for (id in names(all_predictions)) {
  all_predictions.dt <- all_predictions[[id]] %>% # read corresponding data
    right_join(data.frame(month = seq(head(.$month, 1), tail(.$month, 1), 1/freq) ), by = 'month') %>% 
    arrange(desc(-month))
  cat('\n')
  cat('#### ', id, '   \n') 
  
  var.s <- rep(NA_real_, length(all_predictions.dt$r))
  for (t in 1:(length(var.s) - 1)) {
    var.s[t+1] <- var(all_predictions.dt$r[1:t], na.rm = T) # calculate variance estimation
  }
  all_predictions.dt$var.s <- var.s # historical sample variance
  all_predictions.dt <- all_predictions.dt %>%
    left_join(RF, by = "month") %>%
    rename(Rfree = t30ret) %>% # include the risk-free rate
    mutate(Rfree = lag(Rfree)) # choose the lag of that rate
  
  ## obtain the optimal weights on risky portfolios
  weights_stock <- all_predictions.dt
  for (mu in c("r", "hist_mean", "sop_simple", ratio_names)) {
    weights_stock[[mu]] <- (exp(weights_stock[[mu]]) - 1 - weights_stock[["Rfree"]]) / (risk.coef * weights_stock[["var.s"]])
  } ## calculate the optimal weights
  
  cat('\n')
  ## plot the weights
  plot.ts(ts(weights_stock$sop_simple, start = weights_stock$month[1], frequency = freq) %>% window(start = 1986), xlab = NULL, ylab = NULL, main = substitute(paste("The optimal weights on the ", name, " portfolio with ", gamma, " = ", g, sep = ""), list(name = id, g = risk.coef)))
  lines(ts(weights_stock$hist_mean, start = weights_stock$month[1], frequency = freq), col = 2, lty = 3)
  abline(h = c(0, 1, 1.5), col = "greY50", lty = 3)
  legend("bottomright", legend = c("sop_simple", "hist_mean"), col = 1:2, lty = c(1,3), cex = 0.8)
  
  par(mar = c(5, 4, 3, 2))
  plot.ts(ts(weights_stock[c("hist_mean", "sop_simple", ratio_names)],
             start = weights_stock$month[1], frequency = freq) %>% window(start = 1986),
          plot.type = "single", col = 1:length(c("hist_mean", "sop_simple", ratio_names)),
          ylim = c(-1, 3), xlab = NULL, ylab = NULL)
  abline(h = c(0, 1, 1.5), lty = 2, col = "grey50")
  title(main = substitute(paste("The optimal weight on the ", name, " portfolio with ", gamma, " = ", g, sep = ""),
                          list(name = id, g = risk.coef)))
  legend(xpd = T, "bottom", inset = -0.28, lty = 1, lwd = 2, nc = 4, cex = 0.8, bty = "n",
         legend = c("hist_mean", "sop_simpel", ratio_names),
         col = 1:length(c("hist_mean", "sop_simpel", ratio_names)))
  
  ## 8.1 calculate OOS trading performance ----
  ## using upper and lower boundaries
  rp <- all_predictions.dt %>%
    select(month, r, Rfree) %>%
    mutate(r = exp(r) - 1) # change from the log return back to the simple return
  limit.opt = F # the option to set limit or not
  limit.max <- 1.0 # the upper bound of the portfolio weight
  limit.min <- 0 # the lower bound of the portfolio weight
  for (mu in c("hist_mean", "sop_simple", ratio_names)) {
    w <- weights_stock[[mu]]
    if (limit.opt == TRUE) { # whether to pose a upper-lower bound
      w[w > limit.max & !is.na(w)] <- limit.max
      w[w < limit.min & !is.na(w)] <- limit.min
    }
    rp[[mu]] <- w * rp$r + (1 - w) * rp$Rfree # calculate the portfolio returns
  } # 'rp' has the risky+rff trading portfolio returns
  table3.mean <- apply(na.omit(rp[-1]), 2, mean)
  table3.var <- apply(na.omit(rp[-1]), 2, var)
  table3.cer <- table3.mean - 1/2 * risk.coef * table3.var # certainty equivalent returns
  table3.dt <- cbind.data.frame(mean = table3.mean, 
                                variance = table3.var, 
                                CERs = table3.cer,
                                CEGs_annualised = (table3.cer - table3.cer["hist_mean"]) * freq)
  table3.dt
  TABLE3[[id]] <- table3.dt
  
  table3.dt <- round(table3.dt, digits = 6)
  table3.dt$rowname <- rownames(table3.dt)
  table3.dt$portname <- id
  table3.df <- rbind.data.frame(table3.df, table3.dt) # for forming tables
  
  pfm <- na.omit(rp) %>% # simulate its performance
    mutate(buyhold = cumprod(1 + r),
           index.hm = cumprod(1 + hist_mean),
           index.sop = cumprod(1 + sop_simple)) %>% 
    right_join(data.frame(month = seq(rp$month[1], tail(rp$month, 1), 1/freq) ), by = 'month') %>% 
    arrange(desc(-month)) %>% 
    select(month, buyhold, index.hm, index.sop)
  pfm.ts <- ts(pfm[-1], start = pfm$month[1], frequency = freq)
  plot.ts(pfm.ts, plot.type = "single", col = 1:3, log = "y", lty = c(1, 2, 2), xlab = NULL, ylab = "Performance")
  apply(pfm.ts, 2, max)
  # text(x = end(time(pfm.ts)), y = apply(pfm.ts, 2, max), labels = round(apply(pfm.ts, 2, max), digits = 2), cex = 0.5, col = 1:3, pos = 4)
  legend("topleft", legend = colnames(pfm.ts), lty = c(1,2,3), col = 1:3, cex = 0.8)
  if (limit.opt == TRUE) {
    title(main = substitute(paste(id, " (", limit.min, " <= ", omega, " <= ", limit.max, ")", sep = ""),
                            list(id = id, limit.min = limit.min, limit.max = limit.max)))
  } else {
    title(main = paste(id, " (no restrictions)", sep = ""))
  }
  
  ## TABLE-4 Sharpe ratio and Sharpe ratio gains----
  excess_returns <- rp # store the excess returns
  for (mu in c("r", "hist_mean", "sop_simple", ratio_names)) {
    excess_returns[[mu]] <- excess_returns[[mu]] - excess_returns[["Rfree"]]
  } # constructing the excess returns of the portfolio(s)
  names(excess_returns)
  table4.mean <- apply(na.omit(excess_returns)[-1], 2, mean)
  table4.sd <- apply(na.omit(excess_returns)[-1], 2, sd) # variance estimation is the same as in the CEGs. 
  table4.sr <- table4.mean / table4.sd # get the Sharpe ratio
  table4.dt <- cbind.data.frame(mean = table4.mean,
                                sd = table4.sd,
                                SR = table4.sr,
                                SR_annualised = table4.sr * sqrt(freq),
                                SRG_annualised = (table4.sr - table4.sr["hist_mean"]) * sqrt(freq))
  table4.dt
  TABLE4[[id]] <- table4.dt
  
  table4.dt <- round(table4.dt, digits = 6)
  table4.dt$rowname <- rownames(table4.dt)
  table4.dt$portname <- id
  table4.df <- rbind.data.frame(table4.df, table4.dt) 
  
  cat('\n')
}

Market

BH

BL

BM

SH

SL

SM

Table 3 - Certaint equivalent gains

Trading Strategies: certaint equivalent gains

This table shows the out-of-sample portfolio choice results at annual frequencies from predictive regressions and the SOP method. The trading strategy for each portfolio is designed by optimally allocating funds between the risk-free asset and the corresponding risky portfolio. The certainty equivalent return is \(\overline{rp} - \frac{1}{2} \gamma \hat{\sigma}_{rp}^{2}\) with a risk-aversion coefficient \(\gamma = 3\). The annualised certainty equivalent gain (in percentage) is the annual certainty equivalent gain multiplied by the corresponding frequency (e.g. 12 for monthly data).

dt <- table3.df %>%
  filter(rowname %in% c(ratio_names, "sop_simple")) %>%
  select(CEGs_annualised, rowname, portname)

as.data.frame(matrix(dt$CEGs_annualised, byrow = F, nrow = 6, ncol = 7)) %>%
  `colnames<-`(unique(dt$portname)) %>%
  mutate(Variable = unique(dt$rowname)) %>%
  # round(digits = 4) %>%
  as.tbl() %>%
  select(Variable, unique(dt$portname)) %>%
  gt(rowname_col = "Variable") %>%
  tab_header(title = "Table 3 - Trading Strategies: certainty equivalent gains",
             subtitle = paste(str_to_title(freq_name(freq = freq)), " data starts from ", first(data_decompose$month) + k/freq, " and ends in ", last(data_decompose$month), ".", sep = "")) %>%
  fmt_percent(columns = 2:8, decimals = 2)

Table 3 - Trading Strategies: certainty equivalent gains
Annual data starts from Jul 1987 and ends in Jul 2019.
	Market	BH	BL	BM	SH	SL	SM
sop_simple	1.18%	17.55%	−0.36%	4.92%	2.52%	−9.59%	−15.61%
DP	3.25%	9.51%	−1.07%	21.62%	4.62%	−38.49%	−2.52%
PE	0.13%	3.93%	−2.33%	3.70%	3.40%	−2.60%	2.84%
EY	−0.90%	−2.10%	−0.86%	−1.35%	3.52%	9.29%	2.65%
DY	7.00%	8.92%	1.12%	16.36%	5.29%	−4.59%	2.47%
Payout	0.95%	−21.50%	1.20%	−32.57%	1.35%	10.49%	2.22%

Table 4 - Sharpe ratio Gains

Trading Strategies: Sharpe ratio Gains

This table presents the Sharpe ratio of the out-of-sample performance of trading strategies, allocating funds between risk-free and risky assets for each portfolio. The annualised Sharpe ratio is generated by multipling the annual Sharpe ratio by square root of the corresponding frequency (e.g. \(\sqrt{12}\) for monthly data).

dt <- table4.df %>%
  filter(rowname %in% c(ratio_names, "sop_simple")) %>%
  select(SRG_annualised, rowname, portname)

as.data.frame(matrix(dt$SRG_annualised, byrow = F, nrow = 6, ncol = 7)) %>%
  `colnames<-`(unique(dt$portname)) %>%
  mutate(Variable = unique(dt$rowname)) %>%
  # round(digits = 4) %>%
  as.tbl() %>%
  select(Variable, unique(dt$portname)) %>%
  gt(rowname_col = "Variable") %>%
  tab_header(title = "Table 4 - Trading Strategies: Sharpe ratio gains", 
             subtitle = paste(str_to_title(freq_name(freq = freq)), " data starts from ", first(data_decompose$month) + k/freq, " and ends in ", last(data_decompose$month), ".", sep = "")) %>%
  fmt_number(columns = 2:8, decimals = 4) 

Table 4 - Trading Strategies: Sharpe ratio gains
Annual data starts from Jul 1987 and ends in Jul 2019.
	Market	BH	BL	BM	SH	SL	SM
sop_simple	−0.0033	0.2375	−0.0118	−0.0035	−0.1260	−0.1215	−0.0799
DP	−0.0802	−0.0714	−0.0224	0.0620	−0.0104	−0.1921	−0.1139
PE	−0.0878	0.0428	−0.0842	−0.1181	−0.0711	−0.4409	−0.0265
EY	−0.0632	0.0065	−0.0266	−0.0908	−0.0468	−0.1593	−0.0052
DY	0.0869	−0.0326	0.0680	0.0447	0.0077	−0.0642	−0.0347
Payout	0.0134	−0.0553	0.0469	−0.0906	−0.0221	−0.0123	0.0011

Figure 6 - Sensitivity of Certainty Equivalent Gains relative to Risk-Aversion level

This figure presents the out-of-sample portfolio choice results at monthly frequency from bivariate predictive regressions and the SOP method with different levels of risk-aversion. To show that our previous results hold with respect to investors with different levels of risk aversion, we evaluate the changes in certainty equivalent gains with respect to the changes in the level of risk-aversion. The results of the trading strategy reported here are without trading restrictions (as in Table 5), allocating funds between the risk-free asset and the risky equity portfolio. The portfolio choice results are evaluated in the certainty equivalent return with relative risk-aversion coefficient \(\gamma\), with ${\(0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5\)}$. Risky equity portfolios include the market portfolio and six size and book-to-market equity sorted portfolios, BH, BM, BL, SH, SM and SL. The annualised certainty equivalent gain is the monthly certainty equivalent gain multiplied by twelve. The sample period is from February 1966 to December 2019 and the out-of-sample period starts in March 1986.

## 8. calculate the CEGs ----
TABLE3_2 <- data.frame(matrix(nrow = 0, ncol = 0)) # for the CEGs with different risk-averse coefficients.

### 8.1(new) CEGs with different risk-averse coefficients ----
### This part will be recorded in TABLE3_2.
for (id in names(all_predictions)) { # the sensitivity analysis of CEGs wrt. risk-averse coefficient
  print(id); timestamp()
  all_predictions.dt <- all_predictions[[id]] # read corresponding data
  
  var.s <- rep(NA_real_, length(all_predictions.dt$r))
  for (t in 1:(length(var.s) - 1)) {
    var.s[t+1] <- var(all_predictions.dt$r[1:t])
  }
  all_predictions.dt$var.s <- var.s # historical sample variance
  all_predictions.dt <- all_predictions.dt %>%
    left_join(RF, by = "month") %>%
    rename(Rfree = t30ret) %>% # include the risk-free rate
    mutate(Rfree = lag(Rfree)) # choose the lag of that rate
  
  ## record the annualised CEG (certainty equivalent gains)
  table3.ceg_a <- c()
  
  for (risk.coef in seq(from = 0.5, to = 5, by = 0.5)) { # for different level of risk averse
    
    ## obtain the optimal weights on risky portfolios
    weights_stock <- all_predictions.dt
    for (mu in c("r", "hist_mean", "sop_simple")) {
      weights_stock[[mu]] <- (exp(weights_stock[[mu]]) - 1 - weights_stock[["Rfree"]]) / (risk.coef * weights_stock[["var.s"]])
    }
    
    ## calculate OOS trading performance
    ## using upper and lower boundaries
    rp <- all_predictions.dt %>%
      select(month, r, Rfree) %>%
      mutate(r = exp(r) - 1) # convert back to simple returns
    limit.opt = F # the option to set limit or not
    limit.max <- 1.0 # the upper bound of the portfolio weight
    limit.min <- 0 # the lower bound of the portfolio weight
    for (mu in c("hist_mean", "sop_simple")) {
      w <- weights_stock[[mu]]
      if (limit.opt == TRUE) { # whether to pose a upper-lower bound
        w[w > limit.max & !is.na(w)] <- limit.max
        w[w < limit.min & !is.na(w)] <- limit.min
      }
      rp[[mu]] <- w * rp$r + (1 - w) * rp$Rfree # calculate the portfolio returns
    } # 'rp' has the risky-rf trading portfolio returns
    table3.mean <- apply(na.omit(rp[-1]), 2, mean)
    table3.var <- apply(na.omit(rp[-1]), 2, var)
    table3.cer <- (table3.mean - 1/2 * risk.coef * table3.var) # certainty equivalent returns
    table3.ceg <- (table3.cer - table3.cer["hist_mean"]) * freq # annualised certainty equivalent gains
    table3.ceg_a <- c(table3.ceg_a, table3.ceg[['sop_simple']]) # store the annualised CEGs values
  }
  # store the annualised CEGs for the same portfolio with different risk-averse coefficients.
  TABLE3_2 <- rbind.data.frame(TABLE3_2, table3.ceg_a)
}

## [1] "Market"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "BH"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "BL"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "BM"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "SH"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "SL"
## ##------ Wed Mar 27 09:36:25 2024 ------##
## [1] "SM"
## ##------ Wed Mar 27 09:36:25 2024 ------##

colnames(TABLE3_2) <- seq(from = 0.5, to = 5, by = 0.5)
rownames(TABLE3_2) <- names(all_predictions)
write.csv(TABLE3_2, file = "table3_2.csv")

### draw line graphs for these and show the sensitivity.
dt_3.2 <- (TABLE3_2 * 100) %>% # convert to in percentage
  as.matrix() %>%
  c() %>%
  cbind.data.frame(rep(seq(0.5, 5, by = 0.5), each = 7)) %>%
  cbind.data.frame(rep(names(all_predictions), times = 10))
colnames(dt_3.2) <- c("CEGs", "Risk_coef", "Portfolios")
# dt_3.2 # this is the dataset

# jpeg(filename = "table3_2.jpeg", width = 700, height = 500)
ggplot(dt_3.2, aes(x = Risk_coef, y = CEGs, group = Portfolios)) +
  geom_point(shape = 4) +
  geom_line(aes(linetype = Portfolios)) + 
  scale_y_continuous(breaks = seq(-20, 40, by = 10), limits = c(-20, 40)) +
  # scale_y_continuous(breaks = seq(round(min(dt_3.2$CEGs)) - 5, round(max(dt_3.2$CEGs)) + 5, by = 10), limits = c(-40, 20)) +
  ylab(label = "Annualised Certainty Equivalent Gains (in %)") + 
  xlab(label = "Risk-aversion Coefficient") + 
  theme_linedraw()

# dev.off()

Table 5 - MSPE-adjusted Statistic

MSPE-adjusted Statistic

This table presents the MSEP-adjusted Statistics, evaluating the statistical significance of the out-of-sample R-squared statistics of each model in the corresponding portfolio.

See Rapach et al., (2010) and Clark and West (2007) for the detailed procedure.

table5.df <- data.frame()
for (port in names(TABLE5)) {
  pt <- TABLE5[[port]]
  pt$rowname <- rownames(pt)
  pt$portname <- port
  colnames(pt)[4] <- "star"
  table5.df <- rbind.data.frame(table5.df, pt)
}

table5.output <- gt(table5.df, rowname_col = "rowname", groupname_col = "portname") %>%
  fmt_percent(columns = vars(OOS_r.squared, mspe_pvalue), decimals = 2) %>%
  fmt_number(columns = vars(mspe_t), decimals = 4) %>%
  tab_header(title = "Table 5 - MSPE-adjusted Statistic",
             subtitle = paste(str_to_title(freq_name(freq = freq)), " data starts from ", first(data_decompose$month) + k/freq, " and ends in ", last(data_decompose$month), ".", sep = ""))

table5.output

Table 5 - MSPE-adjusted Statistic
Annual data starts from Jul 1987 and ends in Jul 2019.
	OOS_r.squared	mspe_t	mspe_pvalue	star
Market
DP	−27.19%	0.5602	28.97%
PE	−7.11%	1.0960	14.08%
EY	−6.89%	0.7823	22.00%
DY	−4.24%	0.5744	28.49%
Payout	−7.47%	−0.9079	81.45%
SOP	3.83%	1.2086	11.78%
BH
DP	−4.21%	1.1741	12.46%
PE	8.77%	0.7565	22.75%
EY	3.07%	0.5340	29.86%
DY	−2.93%	0.9934	16.41%
Payout	−16.70%	−0.1528	56.02%
SOP	19.58%	1.2851	10.40%
BL
DP	−2.85%	1.0127	15.95%
PE	−2.62%	0.8203	20.92%
EY	0.84%	0.8797	19.29%
DY	5.01%	1.0713	14.62%
Payout	−13.97%	−1.1876	87.80%
SOP	1.84%	0.8678	19.60%
BM
DP	−19.43%	−0.2561	60.02%
PE	−7.39%	1.1765	12.42%
EY	−4.94%	0.8721	19.49%
DY	−8.50%	−0.6969	75.45%
Payout	−4.41%	0.6039	27.51%
SOP	−0.38%	0.5182	30.39%
SH
DP	−13.54%	1.0345	15.45%
PE	−3.61%	0.6932	24.67%
EY	−3.79%	0.2375	40.69%
DY	−4.58%	0.5380	29.72%
Payout	−4.56%	−1.2747	89.40%
SOP	−12.98%	0.5686	28.68%
SL
DP	−86.77%	0.0543	47.85%
PE	−38.25%	0.7416	23.20%
EY	−14.84%	0.4506	32.77%
DY	−13.88%	−0.2351	59.21%
Payout	−10.89%	−0.9767	83.19%
SOP	−16.13%	−0.2284	58.96%
SM
DP	−78.63%	−1.1257	86.55%
PE	−0.38%	0.2756	39.23%
EY	0.39%	0.4473	32.89%
DY	−29.87%	−1.3762	91.07%
Payout	−1.16%	−0.5740	71.50%
SOP	−30.14%	−1.1490	87.05%